%==================================================================================================
% Matlab code to (set) identify the economies of scale and indidivual RICEB measures as discussed in
% "Marital matching, economies of scale and intrahousehold allocations" 
% by Laurens Cherchye, Bram De Rock, Khushboo Surana and Frederic Vermeulen. 
%==================================================================================================


addpath(genpath('C:\Users\u0100562\Downloads\YALMIP-master'));
% replace the directory with the path where YALMIP is downloaded

%% loading the data
clear all
load('psid_cleaned_onlycouples.mat')

%% Bounds on economies of scale and RICEB measures
clc
ub = 0;               % upper bound on degree of publicness of leisure
DivorceCosts =[];     % divorce costs
bounds = [];          % economies of scale
ricebcomp = [];       % RICEB components
a = [];               % A matrix  

while ub <= 0

    temp_id = id;

while isempty(temp_id)==0  
    fprintf('Working on marriage market id %d with Aleisure upper bound as % d\n',temp_id(1), ub);
    %preparing the indices which are in the marriage market being evaluated
    id_ind = find(id==temp_id(1));
    display('Computing Divorce Costs');
    [temp_DivorceCostsM,temp_DivorceCostsF,temp_DivorceCostsMF] = get_divorce_costs(leisurem(id_ind),leisuref(id_ind), hworkm(id_ind),hworkf(id_ind), q(id_ind),wagem(id_ind),wagef(id_ind),nonlabor(id_ind),ub);
    display('Computing bounds');
    [temp_min_scale,temp_max_scale,temp_mina,temp_maxa,temp_min_scalem,temp_min_privatescalem,temp_min_publicscalem,temp_max_scalem,temp_max_privatescalem,temp_max_publicscalem,temp_min_scalef,temp_min_privatescalef,temp_min_publicscalef,temp_max_scalef,temp_max_privatescalef,temp_max_publicscalef] = get_scale(leisurem(id_ind),leisuref(id_ind), hworkm(id_ind), hworkf(id_ind), q(id_ind),wagem(id_ind),wagef(id_ind),nonlabor(id_ind),temp_DivorceCostsM,temp_DivorceCostsF,temp_DivorceCostsMF,ub);
   
    % merging bounds on economies of scale  and RICEB measures
    bounds_temp = zeros(length(id_ind),11);
    bounds_temp(:,1) = temp_min_scale';                        % minimum economies of scale
    bounds_temp(:,2) = temp_max_scale';                        % maximum economies of scale
    bounds_temp(:,3) = temp_max_scale' - temp_min_scale';      % difference in max and min economies of scale
    bounds_temp(:,4) = temp_min_scalem';                       % minimum male RICEB 
    bounds_temp(:,5) = temp_max_scalem';                       % maximum male RICEB
    bounds_temp(:,6) = temp_max_scalem' - temp_min_scalem';    % difference in max and min male RICEB  
    bounds_temp(:,7) = temp_min_scalef';                       % minimum female RICEB
    bounds_temp(:,8) = temp_max_scalef';                       % maximum female RICEB
    bounds_temp(:,9) = temp_max_scalef' - temp_min_scalef';    % difference in max and min female RICEB
    bounds_temp(:,10) = familyid13(id_ind);                    % id of the household 
    bounds_temp(:,11) = ub*ones(length(id_ind),1);             % upper bound on Aleisure used 
    bounds = [bounds; bounds_temp];                         
    
    % merging RICEB components
    ricebcomp_temp = zeros(length(id_ind),11);                  
    ricebcomp_temp(:,1) = temp_min_privatescalem';             % share of private consumption when male RICEB is minimized
    ricebcomp_temp(:,2) = temp_min_publicscalem';              % share of public consumption when male RICEB is minimized
    ricebcomp_temp(:,3) = temp_max_privatescalem';             % share of private consumption when male RICEB is maximized
    ricebcomp_temp(:,4) = temp_max_publicscalem';              % share of public consumption when female RICEB is maximized
    ricebcomp_temp(:,5) = temp_min_privatescalef';             % share of private consumption when female RICEB is minimized
    ricebcomp_temp(:,6) = temp_min_publicscalef';              % share of public consumption when female RICEB is minimized
    ricebcomp_temp(:,7) = temp_max_privatescalef';             % share of private consumption when female RICEB is maximized
    ricebcomp_temp(:,8) = temp_max_publicscalef';              % share of public consumption when female RICEB is maximized
    ricebcomp_temp(:,9) = familyid13(id_ind);                  % id of the household
    ricebcomp_temp(:,10) = ub*ones(length(id_ind),1);          % upper bound on Aleisure used 
    ricebcomp = [ricebcomp; ricebcomp_temp];
    
    % merging the elements of A matrix 
    a_temp = zeros(length(id_ind),12);                     
    a_temp(:,1:5) = temp_mina;                                 % components of A matrix when economies of scale is minimized
    a_temp(:,6:10) = temp_maxa;                                % components of A matrix when economies of scale is maximized
    a_temp(:,11) = familyid13(id_ind);                         % id of the household
    a_temp(:,12) = ub*ones(length(id_ind),1);                  % upper bound on Aleisure used 
    a = [a; a_temp]; 
    
    % merging divorce costs 
    DivorceCosts_minm = min(temp_DivorceCostsMF,[],2);          
    DivorceCosts_minf = min(temp_DivorceCostsMF)';             
    DivorceCosts_avg = (nansum(temp_DivorceCostsMF)' + nansum(temp_DivorceCostsMF,2))./(sum(~isnan(temp_DivorceCostsMF))' + sum(~isnan(temp_DivorceCostsMF),2));
    
    DivorceCosts_temp = zeros(length(id_ind),8);                    
    DivorceCosts_temp(:,1) = DivorceCosts_minm;                  % minimum divorce cost for male's no blocking pairs condition 
    DivorceCosts_temp(:,2) = DivorceCosts_minf;                  % minimum divorce cost for female's no blocking pairs condition  
    DivorceCosts_temp(:,3) = min(DivorceCosts_temp(:,1:2),[],2); % minimum divorce cost for couple's overall no blocking pairs condition 
    DivorceCosts_temp(:,4) = DivorceCosts_avg;                   % average divorce cost for couple's overall no blocking pairs condition
    DivorceCosts_temp(:,5) = temp_DivorceCostsM;                 % divorce cost for male's individual rationality condition 
    DivorceCosts_temp(:,6) = temp_DivorceCostsF;                 % divorce cost for female's individual rationality condition
    DivorceCosts_temp(:,7) = familyid13(id_ind);                 % id of the household
    DivorceCosts_temp(:,8) = ub*ones(length(id_ind),1);          % upper bound on Aleisure used 
    DivorceCosts = [DivorceCosts; DivorceCosts_temp];               
    
    % go to the data of next marriage market
    temp_id(find(temp_id==temp_id(1)))=[];                       % remove these households and go the households in next marriage market
end

fprintf('Computations for ub equal to %d are finished \n', ub);
ub = ub + 0.05;
end


minscale = bounds(:,1);
maxscale = bounds(:,2);
diffscale = bounds(:,3);
minscalem = bounds(:,4);
maxscalem = bounds(:,5);
diffscalem = bounds(:,6);
minscalef = bounds(:,7);
maxscalef = bounds(:,8);
diffscalef = bounds(:,9);

min_aleisurem = a(:,1);
min_aleisuref = a(:,2);
min_ahworkm = a(:,3);
min_ahworkf = a(:,4);
min_aq = a(:,5);
max_aleisurem = a(:,6);
max_aleisuref = a(:,7);
max_ahworkm = a(:,8);
max_ahworkf = a(:,9);
max_aq = a(:,10);

min_privatem = ricebcomp(:,1);
min_publicm = ricebcomp(:,2);
max_privatem = ricebcomp(:,3);
max_publicm = ricebcomp(:,4);
min_privatef = ricebcomp(:,5);
min_publicf = ricebcomp(:,6);
max_privatef = ricebcomp(:,7);
max_publicf = ricebcomp(:,8);

nbm = DivorceCosts(:,1);
nbf = DivorceCosts(:,2);
nbp_min = DivorceCosts(:,3);
nbp_avg = DivorceCosts(:,4);
stabm = DivorceCosts(:,5);
stabf = DivorceCosts(:,6);


save psid_results_onlycouples;




